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ABSTRACT 

We  propose  a  methodology  to  quantify  errors  and  produce  uncertainty  maps  for  satellite-derived  ocean  color  bio-optical 
products  using  ensemble  simulations.  Ensemble  techniques  have  been  used  by  the  environmental  numerical  modeling 
community  to  propagate  initialization,  forcing,  and  algorithm  error  sources  through-out  the  full  simulation  process,  but 
similar  approaches  have  not  yet  been  applied  to  satellite  optical  properties.  Uncertainties  in  retrievals  of  bio-optical 
properties  from  satellite  ocean  color  imagery  are  related  to  a  variety  of  factors,  including  sensor  calibration,  atmospheric 
correction,  and  the  bio-optical  inversion  algorithms.  Errors  propagate,  amplify,  and  intertwine  along  the  processing  path, 
so  it  is  important  to  understand  how'  the  errors  cascade  through  each  step  of  the  analysis,  to  assess  their  impact  and 
identify  the  main  factors  contributing  to  the  uncertainties  in  the  final  products.  Also,  we  are  interested  in  producing 
short-term  forecasts  of  the  bio-optical  property  distributions,  by  coupling  the  satellite  imagery  with  physical  circulation 
models.  So,  in  addition  to  the  uncertainties  in  the  satellite-derived  bio-optical  properties  due  to  the  above-mentioned 
factors,  the  uncertainties  in  the  model  currents  used  to  advect  the  bio-optical  properties  add  another  layer  of  complexity 
to  the  problem.  We  outline  these  processes  and  present  preliminary  results  for  this  approach. 

Keywords:  remote  sensing,  bio-optics,  ocean  color,  error  analysis,  uncertainties,  ensembles,  forecasting 

1.  INTRODUCTION 

Monte-Carlo  methods  have  been  applied  in  meteorology  and  physical  oceanography  to  predict  model  errors  and  to 
improve  weather  and  oceanic  forecasts1,  by  perturbing  the  dominant  sources  of  uncertainties  (e.g.  initial  conditions, 
model  physics,  etc),  but  similar  approaches  have  yet  to  be  applied  to  satellite  optical  properties  and  forecasts.  Typically, 
to  address  uncertainties  in  water  reflectances  and  bio-optical  properties,  satellite  retrievals  are  compared  to  in  situ 
measurements^,  but  this  approach  has  limitations.  There  are  many  errors  and  sources  of  uncertainty  associated  with 
operational  forecasting,  regardless  of  the  forecast  property.  These  include  incomplete  and/or  inaccurate  observations, 
uncertainties  introduced  through  data  assimilation,  and  unresolved  dynamics  and  instabilities.  The  Monte-Carlo 
ensemble  techniques  (multi-analysis  with  a  single  model  or  algorithm)  are  then  based  on  repeated  model  realizations 
with  different  initial  and  boundary  conditions,  atmospheric  forcing,  resolution,  included  data  sets,  bathymetry  grids,  and 
coefficient  values,  for  example.  An  ensemble  reflects  therefore  known  sources  of  forecast  uncertainty,  and  allows  them 
to  be  evaluated'.  There  are  several  way  to  generate  the  parameter  perturbations  for  the  ensembles,  including  random 
bootstrapping  from  normal  distributions  for  surface  drift  problems4,  or  subspace  methods  such  as  the  Ensemble 
Transform  for  initial  field  estimations.  Similar  approaches  can  be  applied  to  physical  or  optical  properties. 

Super-ensemble  techniques  (multi -model  forecasts  or  algorithms)  generalize  the  standard  approach  by  enabling 
comparison  of  different  models5*  .  As  a  result,  each  ensemble  realization  will  no  longer  be  equally  likely  and  bias 
corrections  can  be  obtained  by  looking  to  the  model  skills  at  different  scales.  Using  super-ensembles,  we  can  “blend” 
multiple  algorithms  to  optimize  results.  Statistics  and  metrics  have  been  utilized  (ensemble  mean,  RMS  error,  and 
spread)  to  provide  some  estimate  of  how  well  the  ensembles  capture  “reality,”  thereby  providing  insight  into  the 
underlying  deterministic  processes  and  aiding  decision  support.  Statistical  bias  of  the  models  can  be  determined  from 
multiple  linear  regressions  of  the  forecasts  against  observational  fields.  This  methodology  has  been  successfully  applied 
in  the  ocean  for  surface  drift  problems  and  sound  speed  profile  estimation;  it  can  also  be  used  to  perturb  the  physics  or 
internal  model  parameters  and  allows  fast  tuning  (or  calibration)  of  mesoscale  ocean  nests  to  specific  areas  and  data  sets. 

* Richard.Gould@nrlssc.navy.mil;  phone  1  228  688-5587;  fax  1  228  688-4149;  http://www733l.nrlssc.navy. mil/ 
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Model  and  algorithm  parameters  can  also  be  adjusted  using  ensemble  techniques,  to  optimize  performance  and  forecast 
skill.  Also,  of  particular  interest  to  the  Navy,  ensemble  forecasts  have  been  used  to  provide  probabilistic  forecasts,  which 
in  many  cases  could  be  more  valuable  than  best  guess  forecasts  to  the  warfighter  for  mission  and  contingency  planning. 

There  are  currently  many  new,  advanced  algorithms  in  the  field  of  ocean  color  remote  sensing  (atmospheric  correction, 
bio-optical  inversion,  vertical  profile  models,  advection  schemes  for  optical  forecasts,  electro-optical  sensor  performance 
models)  for  which  the  errors  and  uncertainties  have  not  been  adequately  evaluated  and  characterized.  The  errors 
propagate  through  each  step  in  the  processing,  from  initial  calibration  of  the  satellite  radiance  data  through  the  final 
estimates  of  target  detection/identification.  Estimation  of  uncertainty  in  ocean  color  satellite  imagery  is  needed  for  most 
data  assimilation  schemes.  Melding  of  observations  and  model  predictions  via  data  assimilation  schemes  is  based  on 
their  prior  uncertainty  estimates.  Therefore,  observational  and  model  forecast  error  covariances  are  fundamental 
components  of  most  existing  data  assimilation  approaches  (including  Kalman  filters  and  smoothers,  variational  and 
hybrid  approaches,  etc).  Although  ensemble-based  methodologies  have  been  applied  by  the  modeling  community  to  help 
constrain  model  initial  fields  and  boundary  conditions,  and  to  gain  a  better  understanding  of  uncertainties  associated  with 
model  results,  these  techniques  have  not  been  applied  to  ocean  color  satellite  imagery  . 

Therefore,  a  strong  requirement  exists  to  provide  the  users  of  the  optical  properties,  including  bio-optical  physical 
modelers  who  will  assimilate  the  imagery  and  Navy  mission  planners  and  warfighters,  with  estimates  of  uncertainty  in 
the  derived  satellite  optical  products.  New  measurement  systems  are  now  available  with  optical  instrumentation  (gliders, 
scanfish,  moorings)  and  new  satellite  sensors  are  coming  online  (H1CO,  VI1RS)  to  help  constrain  the  data  ranges  and 
calibrate  the  ensembles,  so  the  time  is  right  to  develop  these  methodologies  for  ocean  color.  Application  of  ensemble 
techniques  to  satellite  ocean  color  processing  will  greatly  improve  the  value  of  the  derived  optical  products  to  the  end- 
users. 

Our  objective  is  to  develop  and  evaluate  methodologies  to  employ  ensemble  techniques  to  satellite  ocean  color  data,  to 
ultimately  improve  atmospheric  correction  and  bio-optical  inversion  algorithms.  By  performing  quantitative  error 
evaluations  and  error  cascading  to  estimate  uncertainties  in  satellite-derived  surface  bio-optical  properties,  this  approach 
will  provide  important  new  information  in  the  form  of  error  maps  and  will  optimize  the  accuracy  and  value  of  the  optical 
properties,  rather  than  providing  only  “best  guess”  optical  products  to  end-users  as  has  been  done  in  the  past.  We  will 
gain  a  better  understanding  of  the  uncertainties  at  all  levels  of  the  processing  of  satellite  ocean  color  imagery.  In 
addition,  the  uncertainty  estimates  will  accompany  optical  fields  assimilated  into  coupled  biophysical  ecological  models, 
improving  downstream  model  forecast  results. 

2.  METHODS 

We  are  developing,  evaluating,  and  applying  ensemble  methodology  to  satellite  ocean  color  imagery.  We  will  extend 
advances  by  the  numerical  modeling  community  to  represent  observational  and  algorithm  error  sources  for  ocean 
ensembles,  and  we  will  evaluate  the  ensemble  representations  and  assess  ensemble  metrics.  The  ensemble  approach  will 
be  applied  at  several  levels:  on  the  initial  satellite  data  (radiances),  on  the  atmospheric  correction  scheme,  on  the  bio- 
optical  inversion  algorithms  used  to  derive  the  optical  properties  (partitioned  absorption  coefficients,  backscattering 
coefficients,  chlorophyll  concentration,  suspended  particulates,  diffuse  attenuation,  and  euphotic  depth),  on  the  second- 
order  optical  products  derived  for  Navy  applications  (diver  visibility,  laser  penetration,  probability  of  target 
detection/identification),  and  finally  on  the  optical  field  forecast  distributions  derived  from  coupling  image  products  with 
circulation  models.  This  will  enable  us  to  examine  the  effects  of  error  propagation  on  the  end-to-end  processing.  Here, 
we  focus  on  uncertainties  due  to  just  two  of  these  aspects,  the  first  and  last  items  above.  First,  we  examine  the  effect  of 
introducing  noise  to  the  top-of-atmosphere  (TO A)  satellite  radiance  values,  and  second,  the  effect  of  perturbations  to  the 
initial  forcing  of  the  physical  circulation  model  used  to  advect  the  bio-optical  properties  (chlorophyll,  in  this  case).  A 
companion  paper  presented  at  this  conference8  addresses  aspects  of  uncertainties  related  to  the  atmospheric  correction; 
future  work  will  address  the  remaining  components.  The  steps  are  illustrated  schematically  in  Figure  1  and  involve 
adapting  methodology  from  numerical  modeling  to  bio-optical  remote  sensing. 

Several  steps  are  required  to  accomplish  these  goals.  We  must  first  calculate  expected  ranges  of  variability  in  satellite 
measurements  of  the  spectral  remote  sensing  reflectances  at  regional  scales  from  an  image  archive.  Also,  existing  data 
sets  collected  during  field  campaigns  will  be  mined  and  assessed  for  suitability  as  training  data  sets.  The  selected  data 
sets  should  contain  as  complete  a  set  of  optical  measurements  as  possible,  covering  all  parameters  estimated  with  the 
satellite  inversion  algorithms.  The  data  sets  should  also  cover  a  wide  range  of  optical  conditions,  from  turbid,  optically- 
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Figure  1.  Schematic  illustrating  the  steps  involved  with  estimating  bio-optical  property  uncertainties  using  an 
ensemble  approach.  Uncertainties  at  each  stage  of  the  processing  contribute  to  the  final,  cumulative  error  map. 


complex  coastal  waters,  to  clear,  open-ocean  waters.  Our  analyses  using  a  wide  variety  and  range  of  in  situ  and  satellite 
optical  data  from  open-ocean  and  coastal  areas  covering  varied  environmental  conditions  will  enable  us  to  form  realistic 
ensembles  that  mimic  real  environmental  variability  and  appropriately  constrain  and  calibrate  the  ensembles.  This  will 
lead  to  accurate  estimates  of  product  uncertainties,  which  are  required  for  assimilation  of  the  bio-optical  products  into 
models.  Initially,  we  will  develop  ensembles  at  a  single  point  to  determine  error  propagation  characteristics.  Then,  we 
will  extend  the  approach  to  an  entire  satellite  image  to  assess  spatial  error  characteristics,  and  to  a  sequence  of  imagery' 
to  examine  temporal  error  characteristics. 

The  satellite  imagery  was  processed  using  the  NRL  Automated  Processing  System  (APS),  which  is  capable  of 
processing  real-time  and  archived  AVHRR,  SeaWiFS,  MODIS,  and  MERIS  satellite  imagery.  APS  is  a  powerful, 
extendable,  image-processing  tool;  it  is  a  complete  end-to-end  system  that  includes  sensor  calibration,  atmospheric 
correction  (with  near-infrared  correction  for  coastal  waters),  and  bio-optical  inversion.  APS  incorporates,  and  is 
consistent  with,  the  latest  NASA  MODIS  code  (SeaDAS)  and  enables  us  to  produce  the  NASA  standard  SeaWiFS  and 
MODIS  products,  as  well  as  Navy-specific  products  using  NRL  algorithms9.  The  automated,  batch  processing  capability 
of  APS  allows  us  to  rapidly  process  a  suite  of  single  point  or  image  ensembles  (based  on  adjusted  radiance  values  and/or 
algorithm  coefficients)  to  facilitate  calculation  of  product  errors. 

To  advect  the  surface  SeaWiFS  satellite  chlorophyll  field  and  produce  24,  48,  and  72-hour  forecast  simulations,  we  used 
currents  derived  from  the  Relocatable  Navy  Coastal  Ocean  Model  (RELO-NCOM).  RELO-NCOM  is  based  on  a 
standardized  development  and  an  efficient  configuration  management  to  facilitate  transitions  of  new  tools  and  real-time 
configurations  of  regional  high  resolution  (order  1  km)  ocean  predictions.  The  physics  and  numerical  procedures  of 
NCOM  are  based  on  the  Princeton  Ocean  Model  (POM)  and  a  Sigma/Z-level  Model  (SZM).  It  solves  a  three- 
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dimensional,  primitive  equation,  baroclinic,  hydrostatic  and  free  surface  system  using  a  cartesian  horizontal  grid,  a 
combination  of  o/z  level  (i.e.,  bottom -following/constant  depth)  vertical  grid  and  implicit  treatment  of  the  free  surface1  A. 
For  mesoscale  real-time  applications,  boundary'  conditions  are  taken  from  an  operational  run  of  the  global  NCOM 
(GNCOM).  The  domain  of  this  particular  experiment  covered  the  entire  Gulf  of  Mexico  (I8N  98  W,  40N  79  W),  from  1 
April  to  30  October  2010.  The  atmospheric  forcing  was  taken  from  the  regional  15km  COAMPS  run  by  the  Fleet 
Numerical  Meteorological  and  Oceanographic  Center  (FNMOC).  Tides  were  introduced  at  the  boundaries  and  through 
local  tidal  potentials.  The  horizontal  grid  spacing  was  set  at  3km  and  used  50  sigma/z  levels  in  the  vertical.  The  model 
assimilates  local  in-situ  observations  along  with  satellite  altimetry  and  sea-surface  temperature  (SST)  data  using  a 
combination  of  model  analysis  and  data;  all  available  observations  from  global  and  local  data  bases  were  assimilated 
over  the  full  period. 

For  the  chlorophyll  forecasts,  a  “pseudo  3-dimensionaP’  Eulerian  advection  scheme  was  used  (without  molecular  or 
turbulent  diffusion  terms).  With  this  approach,  there  are  essentially  two  vertical  layers,  a  1  meter-thick  surface  layer  and 
a  conceptual  deep  layer  to  preserve  continuity  (i.e.,  there  is  vertical  flux  between  the  two  layers,  but  they  move  together 
horizontally).  These  simulations  only  include  current  advection  and  an  assumed  uniform  vertical  chlorophyll  distribution 
based  on  the  surface  values.  Future  enhancements  will  include  addition  of  diffusion  terms,  full  3D  vertical  layering,  and 
the  capability  to  include  more  realistic  vertical  chlorophyll  profiles.  The  forecast  simulations  do  not  include  any 
assimilation  of  in  situ  chlorophyll  data  or  additional  satellite  imagery,  so  currently  the  values  are  unconstrained.  Also, 
with  this  approach,  there  is  an  implicit  assumption  that  the  bio-optical  property'  (chlorophyll)  is  conservative.  Although 
this  is  not  strictly  true,  of  course,  it  may  be  approximately  valid  over  the  short  time  scales  (1-3  days)  that  we  are 
examining,  particularly  in  coastal  areas  where  transport  processes  might  be  expected  to  dominate  biological  processes. 
Therefore,  we  consider  the  optical  properties  to  be  “pseudo-conservative”  tracers  for  our  purposes.  This  allows  us  to 
ignore  growth  and  grazing  terms  and  treat  the  distributional  changes  as  though  they  are  due  entirely  to  dynamical 
processes11. 


3.  RESULTS 

3.1  Top-of-atmosphere  radiance  sensitivity  test 

We  conducted  sensitivity  tests  by  applying  noise  (measurement  uncertainty)  to  the  TOA  satellite  radiance  values.  The 
goal  is  to  add  noise  (random,  or  directional  to  simulate  sensor  drift  or  degradation  over  time)  and  propagate  this  noise 
through  each  step  along  the  processing  path,  to  assess  the  effect  on  uncertainties  in  the  estimated  inherent  optical 
property  (IOP)  retrievals.  First,  we  examine  the  effect  of  adding  noise  to  the  TOA  total  radiance  in  the  748  nm  band 
(L/748)  on  the  calculated  epsilon  value  (e  =  1^748/1^869  for  MOD1S).  Epsilon  is  used  to  select  the  ocean  color 
atmospheric  correction  models  used  to  estimate  the  spectral  aerosol  radiance.  Thus,  errors  in  £  will  lead  to  errors  in 
aerosol  radiance  at  412  nm  (L*412),  and  subsequently  to  errors  in  the  estimates  of  normalized  water-leaving  radiance  at 
412  nm  (nLw412)  and  bio-optical  properties.  Specifically,  we  will  assess  resulting  errors  in  the  absorption  coefficient  at 
412  nm,  a(412),  the  backscattering  coefficient  at  547  nm,  1^547,  and  the  chlorophyll  concentration,  chi.  In  addition, 
these  tests  allow  us  to  assess  how  much  change  in  Lt(748)  is  required  to  change  the  selection  of  the  aerosol  models. 

We  generated  the  test  set  of  TOA  radiance  ensembles  by  drawing  50  samples  from  a  specified  noise  distribution 
described  below.  The  final  results  show  the  variability  (%  difference)  in  the  derived  bio-optical  properties  due  to  the 
addition  of  noise.  We  have  applied  this  approach  at  a  single  pixel;  work  is  underway  to  generate  and  apply  ensembles  at 
each  image  pixel  across  an  entire  ocean  color  image  to  obtain  uncertainty'  maps  for  the  bio-optical  products.  At  this 
point,  we  are  only  considering  uncertainty  due  to  the  addition  of  TOA  noise,  and  no  other  error  sources. 

First,  we  picked  a  random  day,  19  March  2009,  near  Venice,  Italy,  to  use  for  our  single-point  evaluation.  The  selected 
image  pixel  corresponds  to  the  location  of  an  Aerosol  Robotic  NETwork  -  Ocean  Color  (AERONET-OC)  sampling 
platform1",  although  data  from  the  AERONET-OC  are  not  used  in  this  evaluation.  1^(748)  has  a  value  of  0.869713  at  this 
point.  We  then  generated  a  normally-distributed  date  set  of  1^(748)  values  with  the  mean  centered  at  this  value  and  a 
standard  deviation  representing  a  5%  change  (i.e.,  0.043485  radiance  units).  From  this  distribution  we  drew  50  samples 
from  a  total  of  21  bins  spanning  ±3  standard  deviations  (thus  the  bin  width  is  0.013046). 

Thus,  by  applying  noise  to  Lt(748)  we  change  the  value  of  e.  Figure  2  shows  the  resulting  £  values  for  each  of  the  50 
members  in  the  synthetic  Lt(748)  ensemble  data  set.  Any  change  in  Lt(748)  produces  a  linear  change  in  the  calculated  £ 
value.  As  described  above,  the  Lt(748)  samples  were  drawn  from  a  normal  distribution,  so  the  £  distribution  is  normal  as 
well  (depicted  by  the  red  histograms  along  the  X  and  Y  axes  in  Figure  2). 


Proc.  of  SPIE  Vol.  8175  817506-4 


Downloaded  from  SPIE  Digital  Library  on  30  Nov  201 1  to  128  160  1 12  10  Tents  of  Use  http  //spied!  org.’temis 


Figure  2.  Kpsilon  vs.  L,(748)  [red  points,  left  Y  axis],  and  L*(4I2)  vs.  1^748)  [blue  points,  right  Y  axis]  for  the 
test  ensemble  data  set.  Histogram  along  X  axis  illustrates  the  normal  distribution  of  the  1^(748)  data  set; 
histogram  along  the  left  Y  axis  illustrates  the  normal  distribution  of  the  resulting  e  data  set.  The  resulting  L,(412) 
values  (blue  points)  are  not  normally  distributed. 


From  Figure  2  it  is  apparent  that  although  a  linear  change  in  L|(748)  produces  a  linear  change  in  e,  it  does  not  necessarily 
produce  a  linear  change  in  L,(412).  In  this  test,  L„(412)  is  linear  from  approximately  2.2  standard  deviations  below  the 
mean  Lt(748)  value,  to  about  0.2  standard  deviations  above  the  mean.  Above  +0.2  standard  deviations,  L,(412)  remains 
constant  while  e  continues  to  increase.  To  explain  this,  we  examine  the  selection  of  the  bounding  aerosol  models. 

In  the  atmospheric  correction  scheme,  from  a  suite  of  80  aerosol  models,  two  models  are  selected  that  bound  e;  these  two 
sequential  aerosol  models  are  denoted  modmin  and  modmax  (modmin  is  actually  the  higher  number  of  the  two, 
however).  The  aerosol  model  number  indicates  a  specific  relative  humidity  (RH)  and  aerosol  particle  size  fraction  (% 
fine  particle  mode);  the  first  numeral  in  the  model  number  indicates  the  RH  range  and  the  second  numeral  represents  the 
fine  mode  fraction.  Thus,  there  are  ten  size  modes  for  each  RH  range.  For  example,  model  35  represents  a  RH  in  the 
range  75-80%  (indicated  by  the  3)  and  a  fine  mode  fraction  of  20%  (indicated  by  the  5).  The  model  numbering 
conventions,  RH,  and  size  mode  ranges  are  explained  further  in  a  companion  paper  presented  at  this  conference8.  Using 
the  near-infrared  aerosol  radiances  (L*(748,  La(869)),  the  RH,  and  the  fine  mode  fraction,  spectral  La(X)  at  the  sensor 
visible  wavelengths  are  calculated  using  a  radiative  transfer  model11.  Then,  for  a  given  e  obtained  for  each  image  pixel, 
the  spectral  aerosol  radiances  at  the  visible  wavelengths  are  retrieved  from  a  look-up  table,  by  interpolation  between  the 
aerosol  radiances  for  the  two  bounding  models  selected  . 

As  an  example,  for  a  given  e,  suppose  the  bounding  modmin  number  selected  is  35,  then  modmax  should  be  34. 
However,  at  the  top  or  bottom  of  the  model  number  range  for  a  given  RH,  the  modmin  and  modmax  numbers  are  the 
same;  if  a  modmin  of  30  is  selected,  then  modmax  will  also  equal  30,  rather  than  29.  This  is  because  the  aerosol  models 
selected  have  to  stay  within  the  same  relative  humidity'  model  range  (30-39  in  this  example).  The  relative  humidity  is 
selected  before  the  size  fraction,  and  this  prevents  the  models  from  escaping  the  bounds  of  the  current  relative  humidity 
range.  The  same  is  true  at  the  upper  end  of  the  model  index  range. 
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Figure  3A  shows  the  bounding  aerosol  models  selected  when  noise  is  applied  to  U(748)  in  this  sensitivity  test.  For  our 
test  ensemble  data  set,  at  an  e  value  of  approximately  1.19  (roughly  0.2  standard  deviations  above  the  mean  Lt(748) 
value),  modmin  and  modmax  are  both  30,  meaning  that  no  matter  how  much  more  we  increase  Lt(748),  which  in  turn 
produces  a  higher  8  value,  we  cannot  select  new  bounding  aerosol  models.  Because  the  bounding  aerosol  models  are  not 
changing,  we  continue  to  retrieve  the  same  look-up  value  for  La(412),  causing  the  plateau  feature  in  Figure  2.  So,  at  this 
location  on  this  day,  there  is  little  room  for  uncertainty,  such  as  calibration  errors;  variation  greater  than  +0.2  standard 
deviations  for  Lt(748)  will  yield  constant  La(412)  retrievals.  Thus,  in  Figure  2,  as  e  increases,  La(412)  increases  only  to  a 
certain  point;  once  e  reaches  a  certain  point,  it  can  no  longer  be  used  to  compute  two  new  bounding  aerosol  models.  This 
indicates  that  perhaps  modifications  to  the  model  look-up  table  or  radiative  transfer  model  calculations  are  required  to 
improve  the  atmospheric  correction.  We  examine  the  effect  of  an  optimal  model  selection  scheme  in  a  companion  paper 
at  this  conference8. 

The  Ltt(412)  value  subsequently  determines  the  value  of  nLw(412),  with  higher  La  values  leading  to  lower  nU*  values,  so 
the  distribution  of  nLw(412)  in  Figure  3B  mirrors  that  of  La(412)  in  Figure  2.  We  then  propagated  the  50-member 
ensemble  of  Lwn  values  through  the  bio-optical  inversion  algorithms  to  estimate  chlorophyll  (oc3  algorithm14),  a(4 1 2), 
and  ^(547)  (QAA  algorithm15).  The  noise  ensemble  introduced  to  the  Lt(748)  values  led  to  observed  differences  of  -22 
to  +1%,  -75%  to  +8%,  +124%  to  -10%,  for  the  three  bio-optical  properties,  respectively,  from  the  “original”  values 
(those  computed  from  the  image  without  any  noise  applied).  We  are  currently  extending  this  single-point  1^(748) 
ensemble  generation  approach  across  an  entire  ocean  color  image,  to  examine  spatial  variability. 


EPsilon  Lt(748)  (mW/cmJ/MnVsr) 

Figure  3.  Atmospheric  correction  results  for  the  test  ensemble  data  set.  A.  Bounding  model  number  vs.  epsilon.  B. 
nU(412)vs.  U748). 


3.2  Simulated  sensor  degradation 

Ocean  color  sensors  in  space  can  degrade  over  time,  due  to  their  exposure  to  harsh  conditions.  This  leads  to  erroneous 
radiance  measurements  if  left  uncorrected.  For  example,  over  13  years,  the  normalized  radiances  in  the  765  nm  SeaWiFS 
band  showed  a  consistent,  but  non-linear  decrease  of  about  10%,  and  the  865  nm  band  decreased  by  over  20% 16. 
Temporal  calibration  changes  have  been  tracked  for  both  the  Aqua  and  Terra  MOD1S  sensors  as  well17,18.  Although  the 
two  near-infrared  bands  are  used  for  atmospheric  correction,  it  is  also  essential  to  track  the  changes  in  all  the  bands  and 
apply  vicarious  corrections2.  A  systematic  increase  or  decrease  in  the  TOA  radiance  values  recorded  due  to  sensor  drift 
or  degradation  will  lead  to  erroneous  estimates  of  water  bio-optical  properties.  Climate-quality  data  records  (long  term, 
consistent  measurements)  are  required  for  sensor  inter-comparisons  and  to  assess  environmental  trends,  to  ensure  that 
any  observed  data  trend  is  real,  and  not  due  to  sensor  measurement  error.  For  example,  accurate,  long-term  data  records 
are  required  to  determine  whether  global  warming  is  leading  to  any  concomitant  change  in  global  chlorophyll 
concentrations. 
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To  simulate  the  effect  of  sensor  drift  over  time,  we  applied  directional  noise  to  the  1^(765)  values  for  an  example 
SeaWiFS  image.  Figure  4  shows  a  SeaWlFS  chlorophyll  image  of  the  northern  Gulf  of  Mexico  collected  26  April  2010. 
Figure  4A  shows  the  '‘original”  chlorophyll  image  with  no  noise  applied,  and  Figure  4B  shows  the  result  after  applying  a 
1 0.5%  decrease  in  Lt(765)  at  each  pixel.  Application  of  the  noise  led  to  a  decrease  in  the  coastal  chlorophyll  values  and 
an  increase  in  the  offshore  chlorophyll  values.  Keep  in  mind  that  the  chlorophyll  algorithm  is  multi-banded,  so  it  is 
influenced  by  the  spectral  shape  of  the  water-leaving  radiances.  Also,  the  atmospheric  correction  relies  on  the  two  near- 
infra-red  bands,  and  we  are  only  changing  one  of  those  bands.  Therefore,  there  is  not  a  simple  relationship  between  noise 
applied  to  one  atmospheric  correction  band  and  the  resulting  chlorophyll  value. 


Figure  4.  SeaWiFS  chlorophyll  image  covering  the  northern  Gulf  of  Mexico,  26  April  2010.  A.  Standard  product  with  no 
noise  applied.  B.  Result  after  applying  10.5%  decrease  in  Lt(765)  at  each  image  pixel. 


3.3  Chlorophyll  advection  ensembles 

The  chlorophyll  forecast  simulations  were  performed  using  16-member  ensembles  generated  from  the  RELO-NCOM; 
each  ensemble  uses  a  different  forcing  and  initial  condition  (assumed  to  be  the  major  sources  of  error).  The  perturbations 
on  the  initial  conditions  are  set  using  a  technique  called  Ensemble  Transform  that  re-sets  at  each  re-initialization  the 
variances  of  the  perturbations  to  match  our  best  guess  of  the  analysis  error  variances  as  taken  from  the  Navy  Coupled 
Ocean  Data  Assimilation  (NCODA)  system,  after  assimilation  of  observations.  The  atmospheric  forcing  perturbations 
are  taken  from  an  ensemble  population  that  uses  a  space-time  deformation  of  the  deterministic  Coastal 
Ocean/Atmosphere  Mesoscale  Prediction  System  (COAMPS)  runs.  Essentially  this  assumes  COAMPS  has  the  skill  to 
represent  the  major  surface  forcing  features,  but  with  slight  distortions  or/and  delays. 

The  errors  in  the  variable  forecasts  are  determined  by  multiple  sources  of  uncertainty  associated  with  the  model 
initialization  and  boundary  conditions,  numerical  approximations,  modeling  strategies,  impact  of  under-sampling  in  the 
assimilation  process  and  unresolved  scales.  To  address  the  initialization  error  the  best  available  estimate  of  analysis  error 
variance  as  delivered  by  the  data  assimilation  system  is  used  to  transform  forecast  perturbations  into  analysis 
perturbations  by  finding  distinct  linear  combinations  of  forecast  perturbations.  The  analysis  perturbations  are  then  added 
to  the  best  available  analysis  (named  control  run)  to  create  K  initial  states.  These  K  initial  states  are  then  integrated 
forward  in  time  using  the  non-linear  model  to  create  the  next  ensemble  forecast.  This  will  be  then  the  starting  point  for 
the  next  transformation  that  will  generate  the  initial  conditions  for  the  subsequent  ensemble  forecast  once  the  next 
analysis  is  available.  This  technique  does  not  account  for  additional  error  sources  that  could  develop  during  the  forecast 
period  through  the  boundaries  and/or  through  unrepresented  numerical  errors  and  approximations. 

However,  the  NCOM  ensembles  also  represent  the  errors  due  to  uncertain  atmospheric  forcing  by  using  an  ensemble  of 
atmospheric  perturbations  and  allow  each  ocean  run  to  have  an  independent  atmospheric  forcing,  therefore  augmenting 
the  domain  spanned  by  the  ensemble  of  initial  states.  For  this  particular  experiment  the  atmospheric  ensemble  was 
generated  by  atmospheric  perturbations  using  spatial ly-varying  time  shifts  of  the  forecast,  with  a  choice  of  parameters  to 
provide  a  well  developed  spread  of  atmospheric  perturbations.  This  method  is  mostly  adequate  when  predicted 
atmospheric  fields  contain  the  forecast  features  of  interest,  but  they  are  misplaced  in  space  and  time. 

These  ensembles  resulting  from  combining  perturbed  initialization  and  atmospheric  forcing  are  used  to  predict  how 
uncertainties  of  the  ocean  fields  will  evolve  in  space  and  time.  Typical  implementations  of  ensemble  runs  use  order  40 
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independent  members  but  this  particular  experiment  only  had  16  runs.  The  ensembles  thus  represent  a  suite  of  model 
current  vectors  with  varying  initial  and  forcing  conditions  that  simulates  uncertainty  in  the  model  forecasts.  The  16 
model  current  ensembles  are  then  used  to  create  24,  48,  and  72  hour  forecasts  of  the  surface  chlorophyll  distribution  by 
advecting  the  satellite  chlorophyll  image.  Thus,  16  chlorophyll  forecasts  are  generated  from  the  16  model  ensemble  runs, 
along  w  ith  a  control  run  (no  perturbation  of  the  model  initial  forcing  conditions). 


Figure  5  shows  the  results  for  the  48  hour  forecast  (24  and  72-hour  results  not  shown).  A  cloud-free  SeaWiFS  image 
collected  on  26  April  2010  covering  the  Mississippi  Bight  in  the  northern  Gulf  of  Mexico  was  selected  for  the  test  (Fig. 
5A).  The  forecast  from  the  control  run  (Fig.  5B)  shows  distinct  chlorophyll  features  (due  to  the  fact  that  we  are  not 
presently  accounting  for  diffusion  terms  in  the  model);  these  features  are  less  distinct  (more  “blurred”  or  smoothed)  in 
the  forecast  ensemble  mean  (Fig.  5C).  The  normalized  standard  deviation  of  the  chlorophyll  forecast  (standard  deviation 
of  the  16  ensemble  chlorophyll  values  at  each  pixel  divided  by  the  ensemble  mean)  is  shown  in  Figure  5D.  Only  large- 
scale,  offshore  currents  are  constrained  in  the  model  ensembles  (through  assimilation  of  satellite  SSH  into  RELO- 
NCOM);  coastal  currents  are  not  constrained  (i.e.,  no  CODAR  currents  assimilated).  This  contributes  to  the  high 
variability  in  the  chlorophyll  forecasts  observed  in  the  frontal  zone  area  west  of  the  Mississippi  River  delta,  indicating 
the  region  where  the  model  currents  exhibit  the  largest  variations.  Without  constraining  model  current  velocities,  the 
uncertainty  in  the  model  predictions  of  the  chlorophyll  tracer  is  correspondingly  high. 


Figure  5.  48-hour  chlorophyll  forecast  simulation,  northern  Gulf  of  Mexico.  A.  Initial  SeaWiFS  scene  from  26 
April  2010.  B.  RELO-NCOM  forecast  for  the  control  run  (no  perturbation  of  initial  conditions).  C.  Mean  forecast 
from  the  1 6-member  ensemble  runs.  D.  Normalized  standard  deviation  of  the  ensemble  forecasts. 


4.  CONCLUSIONS 

We  have  outlined  an  approach  to  apply  ensemble  techniques  to  estimate  uncertainties  in  satellite-derived  bio-optical 
properties.  We  addressed  uncertainties  in  the  TOA  radiances  recorded  at  the  sensor,  and  their  propagation  to  the  final 
bio-optical  product  estimates.  Based  on  a  limited  set  of  single-point  error  characteristics,  initial  results  indicate  that 
introduction  of  noise  to  the  TOA  radiances  at  the  level  of  approximately  -10%  to  +10%  leads  to  corresponding  errors  of 
-75%  to  +8%  in  the  derived  a(412)  values,  +124%  to  -10%  in  the  derived  bb(547)  values,  and  -22  to  +1%  in  the  derived 
chlorophyll  values.  Initial  work  to  extend  this  single-point  ensemble  generation  approach  spatially  across  an  ocean  color 
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image  is  underway.  We  advected  a  satellite  chlorophyll  field  using  a  suite  of  16  current  ensembles  derived  from  a 
circulation  model,  and  produced  mean  and  standard  deviation  fields  corresponding  to  a  24,  48,  and  72-hour  forecasts. 
These  results  demonstrate  the  utility  of  ensemble  techniques  to  estimate  uncertainties  of  bio-optical  properties  derived 
from  satellite  ocean  color  imagery. 

5.  FUTURE  WORK 

We  will  continue  to  employ  ensemble  techniques  to  address  uncertainties  and  their  propagation  at  each  stage  along  the 
satellite  ocean  color  processing  path.  We  will  continue  to  apply  and  evaluate  the  single-point  ensemble  development 
approach  spatially  across  an  entire  image,  to  obtain  spatial  uncertainty  maps  for  the  derived  bio-optical  properties.  In 
addition,  we  will  examine  natural  bio-optical  variability,  using  in  situ  measurements  and  archived  ocean  color  imagery, 
to  ensure  that  the  ensembles  we  develop  represent  reality.  We  will  compare  the  ensemble  probability  distribution 
functions  (pdfs)  to  those  from  the  measured  in  situ  values  and  adjust  the  ensembles  through  application  of  ensemble 
transform  techniques,  and  we  will  apply  metrics  to  evaluate  ensemble  skill.  We  will  begin  to  explore  the  effect  of  the 
bio-optical  algorithm  selection  on  the  satellite  bio-optical  products,  and  we  will  apply  ensemble  techniques  to  optimize 
the  bio-optical  inversion  algorithm  coefficients.  For  example,  several  of  the  bio-optical  algorithms  include  coefficients 
that  are  applied  globally  but  were  estimated  empirically  using  limited  data  sets  (for  example,  the  slope  of  CDOM 
spectral  absorption  and  the  spectral  backscattering  coefficient).  These  coefficients  can  be  adjusted  on  a  regional  basis 
using  local  in  situ  training  data  sets  in  ensembles.  Finally,  we  will  employ  super-ensemble  techniques  to  “blend”  several 
bio-optical  algorithms  to  achieve  optimal  regional  performance  (lower  errors).  For  example,  there  are  multiple 
algorithms  available  for  bio-optical  inversion,  including  the  QAA,  GSM1,  PML,  and  others,  and  we  will  examine  these 
algorithms  to  assess  performance  in  different  water  types  (open-ocean,  coastal,  turbid).  For  the  bio-optical  forecasts,  we 
will  attempt  to  assimilate  additional  chlorophyll  imagery  (eg.,  satellite  imagery  from  the  days  following  the  initial 
image),  and  we  will  assess  forecast  skill.  Finally,  we  will  produce  uncertainty  maps  for  the  satellite-derived  bio-optical 
properties  and  their  short-term  forecasts,  by  combining  the  errors  propagated  along  the  entire  processing  path. 
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